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Abstract.  Eddy  current  models  have  matured  to  such  a  degree  that  it  is  now 
possible  to  simulate  realistic  nondestructive  inspection  (NDI)  scenarios.  Models 
have  been  used  in  the  design  and  analysis  of  NDI  systems  and  to  a  limited  extent, 
model-based  inverse  methods  for  Nondestructive  Evaluation  (NDE).  The  science 
base  is  also  being  established  to  quantify  the  reliability  of  systems  via  Model- 
Assisted  Probability  of  Detection  (MAPOD).  In  realistic  situations,  it  is  more 
accurate  to  treat  the  input  model  variables  as  random  variables  rather  than 
deterministic  quantities.  Typically  a  Monte-Carlo  simulation  is  conducted  to 
predict  the  output  of  a  model  when  the  inputs  are  random  variables.  This  is  a 
reasonable  approach  as  long  as  computational  time  is  not  too  long;  however,  in 
most  applications,  introducing  a  flaw  into  the  model  results  in  extensive 
computational  time  ranging  from  hours  to  days,  prohibiting  Monte-Carlo 
simulations.  Even  methods  such  as  Latin-Hypercube  sampling  do  not  reduce  the 
number  of  simulations  enough  for  reasonable  use.  This  paper  presents  the 
Probabilistic  Collocation  Method  as  a  non-intrusive  alternative  to  other 
uncertainty  propagation  techniques. 
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Introduction 

Recent  model  development  efforts  have  demonstrated  significant  advances  in  the 
simulation  of  eddy  current  inspection  scenarios  [1-3].  This  work  is  a  critical 
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component  of  a  strategy  to  use  modeling  and  simulation  in  NDE  to  reduce  the  cost  and 
experimental  burden  of  Probability  of  Detection  (POD)  studies,  characterize  damage, 
and  optimize  probe  design  and  inspection  configuration. 

The  authors  have  investigated  the  utility  of  several  numerical  methods  for 
electromagnetic  NDE.  These  methods  include  analytical  methods,  the  finite  difference 
method,  the  finite  element  method  (FEM),  meshless  FEM,  the  boundary  element 
method  (BEM),  and  the  volume  integral  method.  In  addition  to  these  methods,  hybrid 
methods  are  also  being  developed  [2].  While  incremental  advances  in  computational 
efficiency  are  expected,  the  future  development  of  model-assisted  probability  of 
detection  (MAPOD)  protocols  and  model-based  inversion  schemes  should  not  depend 
on  a  revolution  in  computational  efficiency.  One  technical  capability  necessary  for  the 
realization  of  MAPOD  and  deployable  model-based  inverse  methods  is  the  propagation 
of  uncertainty  in  input  variables.  In  this  work,  variable  refers  to  a  measurable  quantity, 
and  parameter  refers  to  a  quantity  that  must  be  estimated  such  as  liftoff. 

Stochastic  numerical  methods  have  been  developed  to  propagate  uncertainty 
efficiently  without  relying  solely  on  advances  in  computational  power.  Significant 
advances  in  stochastic  computation  have  been  made  especially  in  the  last  20  years. 
This  paper  briefly  traces  the  development  of  stochastic  numerical  methods  pertinent  to 
NDE  modeling  and  simulation.  The  probabilistic  collocation  method  (PCM)  is 
introduced  and  applied  to  two  eddy  current  forward  models. 


1.  Probabilistic  Collocation  Method 

Spectral  representation  of  uncertainty  involves  decomposing  a  random  variable  into 
deterministic  and  stochastic  components.  Following  the  work  of  Norbert  Wiener  on 
Homogeneous  Chaos  [4],  Cameron  and  Martin  pointed  out  that  any  second-order 
functional  of  Brownian  motion  can  be  expressed  as  a  mean-square  convergent  series  in 
terms  of  infinite-dimensional  Hermite  polynomials  in  Gaussian  variables  [5].  Any 
random  process  with  finite  second-order  moments,  (which  encompasses  most  physical 
phenomena)  can  be  represented  by  a  Polynomial  Chaos  expansion  using  Hermite 
polynomials  and  will  converge  according  to  the  Cameron-Martin  theorem.  If  this 
random  process  is  Gaussian,  the  convergence  is  exponential.  An  important  property  of 
Hermite  polynomials  is  that  they  are  orthogonal  with  respect  to  Gaussian  probability 
measure.  In  fact  for  Gaussian  processes,  optimal  convergence  is  achieved  with  the  use 
of  Hermite  polynomials  because  the  weighting  function  is  the  same  as  the  probability 
density  function  (PDF)  for  the  Gaussian  distribution.  To  illustrate,  the  definition  of  an 
orthogonal  polynomial  with  respect  to  a  given  weighting  function  w(x)  is  given  by: 

b 

(/  •  g)„  =  j  f(x)g{x)w{x)dx  =  0  (1) 


The  inner  product  of  functions  f(x)  and  g(x)  is  zero  if  they  are  orthogonal  with  w(x). 
The  only  polynomial  functions  that  satisfy  this  for  the  case  where  the  weighting 
function  is  equal  to  the  PDF  of  a  standard  normal  distribution  are  Hermite  polynomials. 


Ghanem  and  Spanos  combined  this  concept  with  the  finite  element  method  to 
simulate  problems  in  solid  mechanics  [6].  For  non-Gaussian  processes,  the 
convergence  is  not  optimal,  but  fortunately  there  are  other  orthogonal  polynomial  sets 
for  other  general  types  of  random  variables  that  converge  exponentially  [7]. 

In  work  mentioned  thus  far,  the  formulations  are  often  intrusive,  meaning  that 
alternation  of  deterministic  code  is  necessary  for  implementation.  A  non-intrusive 
method  was  introduced  in  the  mid- 1990’s  [8].  Essentially,  a  reduced  order  polynomial 
form  of  the  model  is  derived  with  minimal  simulations.  The  input  values  for  these 
simulations  are  derived  from  the  roots  of  orthogonal  polynomials  which  are  dependent 
only  on  the  input  distributions.  Hermite  polynomials  are  associated  with  Gaussian 
distributions  as  described  in  the  preceding  paragraph.  The  weighting  function  for  the 
Uniform  distribution  is  simply  a  constant  on  an  interval  bounded  by  [a,b].  Legendre 
polynomials  are  orthogonal  with  this  weighting  function. 


2.  Case  Study  1. 

Application  of  the  Probabilistic  Collocation  Method  is  model  independent,  but  the 
decision  of  what  order  and  whether  to  include  interaction  terms  does  depend  on  the 
complexity  of  the  mapping  of  input  variables  to  the  output  response.  The  orthogonal 
polynomials  and  their  associated  roots  are  dependent  on  the  type  of  input  distribution 
and  its  parameters.  For  this  case  study,  the  teaml5B  workshop  problem  presented  by 
Steve  Burke  is  used  [9].  All  simulations  are  conducted  using  VIC-3D®.  The  only 
variable  that  is  changed  in  this  example  is  the  frequency,  which  is  500  Hz  here  instead 
of  7  kHz.  In  this  problem,  the  probe  is  scanned  along  the  notch  and  the  real  and 
imaginary  components  of  the  impedance  are  measured.  A  peak  in  the  magnitude  of  the 
impedance  is  observed  when  the  center  of  the  coil  translated  along  the  notch  at  a 
distance  of  -12  mm  and  12  mm  away  from  the  origin  which  is  the  center  of  the  notch. 

The  objective  of  this  exercise  is  to  predict  the  probability  density  function  (PDF) 
of  the  peaks  of  the  real  and  imaginary  components  of  the  impedance  when  the  liftoff  of 
the  coil  is  uniformly  distributed  between  1.53  mm  and  2.53  mm,  and  the  depth  of  the 
notch  is  normally  distributed  with  a  mean  of  5  mm  and  standard  deviation  of  1  mm.  In 
the  original  problem  the  coil  liftoff  is  2.03  mm  and  the  depth  of  the  flaw  is  5  mm. 

The  first  step  is  to  derive  orthogonal  polynomials  with  respect  to  the  input 
distributions.  From  now  on  the  parameter  liftoff  will  be  denoted  by  ‘A’,  and  notch 
depth  will  be  denoted  by  ‘B’.  For  liftoff,  Legendre  polynomials  associated  with  that 
particular  distribution  and  parameter  set  are  shown  in  equation  2.  For  the  depth,  the 
Hermite  polynomials  are  shown  in  equation  3. 

L{(A)  =  A-  2.03 
L2{A)  =  A2  —4.06/1  +  4.0376 
L2(A)=  A3 -6.09  A2  +  12.2127/1-8.0609 
L4(A)  =  A4  -8.1243  +  24.51 1 1 A2  —32.5917^4  +  16.1041 
L5(A)  =  A5  - 10. 1544  +  40.93 12 A3  -81 .9626 A2  +  8 1 .48994 - 32.1796 


(2) 


h2(C)=C2~  i 
h3(C)=C3-  X 


//4(o=r-6^2+ 3 

H5(£)  =  £5  -10i^3  + 15£" 


(3) 


Figure  1 .  Change  in  the  (a)  real  and  (b)  imaginary  part  of  the  impedance  due  to  the 
plate  and  flaw. 


Note  that  since  the  mean  of  the  Gaussian  distribution  is  5,  parameter  B  is  translated  by 
5. 


B  =  5  +  Hl(0  (4) 

Next,  a  first  order  model  might  be  used  to  approximate  the  real  model.  In  this 
case,  there  are  three  unknowns  X0,  X1?  and  X2  as  shown  in  equation  5.  Since  we  are 
interested  in  both  the  resistance  and  reactance,  this  process  is  done  separately  for  each. 
Three  simulations  are  required  to  solve  for  the  three  coefficients.  Four  pairs  of 
collocation  points  for  the  inputs  of  the  simulations  are  available  from  the  roots  of  the 
2nd  order  polynomials.  Three  pairs  are  selected  to  solve  for  the  three  unknowns. 
and  Hi  are  simply  evaluated  at  the  three  pairs  of  collocation  points.  This  is  all  that  is 
needed  for  the  a  1st  order  approximate  model.  Then  a  Monte-Carlo  simulation  can  be 
conducted  on  this  simple  model  to  calculate  the  PDF  of  the  resistance  and  reactance. 
To  check  the  accuracy  of  this  approximation,  the  next  higher  order  collocation  points 
are  necessary.  The  roots  of  the  3rd  order  orthogonal  polynomials  are  used  and  six 
simulations  using  the  full  model  are  conducted  and  compared  with  the  approximate 
model.  Ultimately  nine  simulations  are  necessary  for  construction  and  evaluation  of 
the  1st  order  approximate  model.  Fortunately  the  same  six  simulations  used  for  the 
error  check  of  the  1st  order  model  can  be  used  to  construct  the  2nd  order  approximate 
model.  The  4th  order  orthogonal  polynomials  are  then  used  to  check  the  accuracy  of  the 
2nd  order  approximate  model  requiring  eight  more  simulations  for  a  total  of  fourteen 
simulations.  The  authors  have  generally  found  that  a  3rd  order  approximate  model  with 
one  interaction  term  is  usually  sufficient  to  predict  output  PDF’s  with  good  accuracy 
for  forward  eddy  current  simulations,  but  there  may  be  situations  where  higher  order 
approximations  with  more  interaction  terms  are  needed. 


Y  =  X0+  XxLx  (A)  +  X2Hx  ( B ) 


(5) 


The  form  of  the  3rd  order  reduced  model  as  described  in  equation  6  has  eight 
coefficients,  one  of  which  is  an  interaction  term. 


Y  =  X0  +  XlLl(A)+  X2Hy($)+  X3L2(A)+  X4H2(£) 
+  X5L3(A)+  X6H3(f)+X7L1(A)H1(f) 


(6) 


To  solve  for  the  eight  unknowns,  eight  simulations  must  be  run.  The  roots  of  the 
4th  order  Legendre  and  Hermite  polynomials  are  used  to  select  the  input  values  for  A 
and  B  in  these  simulations.  This  is  done  separately  for  the  real  and  imaginary 
components  of  the  impedance.  Once  again,  everything  is  known  except  for  the  X 
coefficients.  The  resulting  system  of  equations  is  then  solved  and  the  coefficients  are 
listed  in  Table  1  for  both  the  real  and  imaginary  components  of  the  impedance. 
Incidentally,  X0  is  the  expected  value,  which  is  one  of  the  primary  benefits  of  using 
orthogonal  polynomials.  As  expected,  the  interaction  term  X7  is  close  to  zero  in  both 
cases. 

To  check  the  error  of  this  third  order  approximation,  the  roots  of  the  5th  order 
polynomials  are  used  to  select  collocation  points.  Ten  more  simulations  are  needed  for 
comparison  with  the  reduced  order  model  predictions.  The  results  of  the  simulations 
and  the  results  of  the  reduced  model  are  displayed  in  Table  2. 


Table  1.  Coefficients  for  3rd  order  approximation  model. 


Model  term 

Real 

Imaginary 

X0 

0.9693 

-2.4132 

Xi 

-0.2209 

0.4386 

X2 

0.0012 

0.0072 

X3 

0.0284 

-0.0479 

X4 

0.0004 

-0.0011 

X5 

-0.0050 

0.0077 

X6 

-0.0001 

0.0003 

X7 

-0.0003 

-0.0018 

Table  2.  Evaluations  of  original  and  approximation  model  for  error  calculation. 


Y  real 

Y  real 

Y  imaginary 

Y 

1  imaginary 

1.027 

1.0270 

-2.5444 

-2.5415 

1.0284 

1.0280 

-2.5272 

-2.5295 

1.0306 

1.0305 

-2.5202 

-2.5214 

0.9656 

0.9656 

-2.4219 

-2.4195 

0.9670 

0.9666 

-2.4059 

-2.4081 

0.9690 

0.9690 

-2.3994 

-2.4007 

0.9084 

0.9084 

-2.3066 

-2.3043 

0.9098 

0.9092 

-2.2918 

-2.2937 

0.9117 

0.9116 

-2.2857 

-2.2869 

0.8774 

0.8778 

-2.2073 

-2.2053 

The  residual  is  defined  as  the  difference  between  the  full  and  reduced  model  in 
equation  7,  and  the  sum  of  squares  of  the  residuals  is  defined  in  equation  8.  The 
relative  sum  of  squares  of  the  residuals  is  used  to  quantify  the  error  in  the 


approximation.  The  relative  sum  of  squares  of  the  residuals  is  calculated  by  simply 
dividing  by  the  expected  values  as  shown  in  equation  9.  These  quantities  are  summed 
over  the  model  results  for  each  collocation  point. 


s  =  Y-Y 


(7) 


ssr  = 


(8) 


E[Y] 

The  errors  for  the  real  and  imaginary  parts  are  given  by:  ssr( resistance)  = 
2.9428xl0'4,  mr(resistance)  =  3.0359xl0'4  ^r(reactance)  =  0.0020,  and 
mr(reactance)  =  8.4790xl0~4. 

Now  there  are  two  3rd  order  polynomials  with  good  accuracy  for  the  real  and 
imaginary  components  of  the  impedance.  Monte  Carlo  simulations  can  now  be  run 
with  ease  on  this  3rd  order  approximate  model.  The  PDF  for  the  real  and  imaginary 
component  of  the  impedance  is  shown  in  Figure  2(a)  and  2(b)  respectively. 


Figure  2.  PDF  of  the  (a)  real  and  (b)  imaginary  component  of  the  impedance. 


3.  Case  Study  2. 

A  second  case  study  is  presented  for  the  characterization  of  a  gap  between  multiple 
plates,  representing  the  problem  of  characterizing  corrosion  at  the  faying  surface  of  a 
multilayer  aircraft  structure.  A  diagram  of  the  problem  is  shown  in  Figure  3  including 
material  properties  and  dimensions.  In  this  hypothetical  study,  the  thickness  and 


conductivity  of  the  layer  are  treated  as  random  variables.  The  conductivity  of  the  top 
plate  is  assumed  to  be  a  normal  random  variable  with  a  mean  value  of  1.876><107  S/m 
and  standard  deviation  of  3. 3  84x1 06  S/m  as  shown  in  Figure  4(a).  The  liftoff  is  defined 
as  a  uniform  random  variable  with  the  range  of  0.55  +/-  0.05  mm  as  shown  in  Figure 
2(b).  VIC-3D®  was  used  to  generate  the  simulated  results  for  the  study.  The  gap 
between  the  plates  was  fixed  at  0.3  mm  and  the  frequency  for  the  study  was  set  to  4.0 
kHz.  Both  1st  order  and  2nd  order  PCM  models  were  evaluated  for  this  case  requiring 
only  3  and  6  model  calls  respectively.  Monte  Carlo  simulations  were  also  run  with 
100,000  model  calls  used  in  the  study. 


dz 

0.81  mm 

dh 

6.55  mm 


1 . . 


.  . 1 


a!  =  1.876  E  107  S/m 
aair  =  0.0  107  S/m 
g2  =  1.856  E  107  S/m 


Figure  3.  Diagram  of  case  study  2  problem. 


Figure  4.  (a)  Conductivity  and  (b)  liftoff  input  distributions  with  ouput  (c)  resistance 
and  (d)  reactance  output  ID  distributions  and  (e)  joint  distribution  calculated  using 
PCM.  (f)  A  joint  distribution  calculated  using  Monte  Carlo  simulations  is  provided  for 
comparison. 


Results  for  the  2nd  order  PCM  model  are  presented  in  Figures  4(c)-(f).  The  results 
in  Figures  4(c)-(d)  show  that  the  resistance  is  basically  distributed  uniformly  and  the 
reactance  is  basically  distributed  as  a  normal  random  variable.  This  relationship  is 
expected  based  on  the  different  effect  of  liftoff  and  conductivity  changes  in  the 
impedance  plane.  A  comparison  of  the  joint  distribution  calculated  using  both  PCM 
and  Monte  Carlo  methods  are  shown  in  Figures  2(e)  and  2(f)  respectively.  By 
observation,  good  agreement  was  achieved  using  PCM  with  Monte  Carlo  methods 


while  only  requiring  6  model  calls  here.  Values  for  the  error  between  the  PCM  and 
Monte  Carlo  simulations  for  the  1st  order  PCM  model,  ssr  and  rssr,  were  43.58  and 
0.3177  respectively.  Error  values  for  the  2nd  order  PCM  model,  ssr  and  rssr,  were 
15.49  and  0.1129  respectively.  Thus,  the  2nd  order  model  provides  a  64.5%  reduction 
in  the  error  between  the  PCM  and  Monte  Carlo  simulations  with  respect  to  the  1st  order 
model.  In  conclusion,  the  PCM  approach  is  highly  efficient  at  propagating  parameter 
variation  through  models  with  minimal  model  calls. 


4.  Summary 

A  non-intrusive  approach  to  propagating  random  inputs  in  eddy  current  forward 
models  was  presented.  The  selection  of  collocation  points  in  this  method  depends  on 
the  type  and  parameters  of  the  input  statistical  distribution  and  not  the  model  itself. 
Prediction  of  the  PDF  can  be  achieved  with  fewer  simulations  than  traditional 
uncertainty  propagation  methods.  Future  work  will  investigate  using  the  method  with 
more  input  variables  and  also  for  inverse  problems.  The  robustness  of  this  method  for 
inverse  problems  will  likely  be  much  more  sensitive  to  choice  of  order  and  interaction 
terms. 
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